home *** CD-ROM | disk | FTP | other *** search
/ io Programmo 60 / IOPROG_60.ISO / soft / c++ / gsl-1.1.1-setup.exe / {app} / src / interpolation / test.c < prev   
Encoding:
C/C++ Source or Header  |  2002-04-18  |  7.1 KB  |  252 lines
  1. /* interpolation/test.c
  2.  * 
  3.  * Copyright (C) 1996, 1997, 1998, 1999, 2000 Gerard Jungman
  4.  * 
  5.  * This program is free software; you can redistribute it and/or modify
  6.  * it under the terms of the GNU General Public License as published by
  7.  * the Free Software Foundation; either version 2 of the License, or (at
  8.  * your option) any later version.
  9.  * 
  10.  * This program is distributed in the hope that it will be useful, but
  11.  * WITHOUT ANY WARRANTY; without even the implied warranty of
  12.  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  13.  * General Public License for more details.
  14.  * 
  15.  * You should have received a copy of the GNU General Public License
  16.  * along with this program; if not, write to the Free Software
  17.  * Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
  18.  */
  19.  
  20.  
  21. /* Author:  G. Jungman
  22.  */
  23. #include <config.h>
  24. #include <stddef.h>
  25. #include <stdlib.h>
  26. #include <stdio.h>
  27. #include <math.h>
  28. #include <gsl/gsl_test.h>
  29. #include <gsl/gsl_errno.h>
  30. #include <gsl/gsl_interp.h>
  31.  
  32.  
  33. int
  34. test_bsearch(void)
  35. {
  36.   double x_array[5] = { 0.0, 1.0, 2.0, 3.0, 4.0 };
  37.   size_t index_result;
  38.   int status = 0;
  39.   int s;
  40.  
  41.   /* check an interior point */
  42.   index_result = gsl_interp_bsearch(x_array, 1.5, 0, 4);
  43.   s = (index_result != 1);
  44.   status += s;
  45.   gsl_test (s, "simple bsearch");
  46.  
  47.   /* check that we get the last interval if x == last value */
  48.   index_result = gsl_interp_bsearch(x_array, 4.0, 0, 4);
  49.   s = (index_result != 3);
  50.   status += s;
  51.   gsl_test (s, "upper endpoint bsearch");
  52.  
  53.   /* check that we get the first interval if x == first value */
  54.   index_result = gsl_interp_bsearch(x_array, 0.0, 0, 4);
  55.   s = (index_result != 0);
  56.   status += s;
  57.   gsl_test (s, "lower endpoint bsearch");
  58.  
  59.   /* check that we get correct interior boundary behaviour */
  60.   index_result = gsl_interp_bsearch(x_array, 2.0, 0, 4);
  61.   s = (index_result != 2);
  62.   status += s;
  63.   gsl_test (s, "degenerate bsearch");
  64.  
  65.   /* check out of bounds above */
  66.   index_result = gsl_interp_bsearch(x_array, 10.0, 0, 4);
  67.   s = (index_result != 3);
  68.   status += s;
  69.   gsl_test (s, "out of bounds bsearch +");
  70.  
  71.   /* check out of bounds below */
  72.   index_result = gsl_interp_bsearch(x_array, -10.0, 0, 4);
  73.   s = (index_result != 0);
  74.   status += s;
  75.   gsl_test (s, "out of bounds bsearch -");
  76.  
  77.   return status;
  78. }
  79.  
  80.  
  81.  
  82. typedef double TEST_FUNC (double);
  83. typedef struct _xy_table xy_table;
  84.  
  85. struct _xy_table
  86.   {
  87.     double * x;
  88.     double * y;
  89.     size_t n;
  90.   };
  91.  
  92. xy_table make_xy_table (double x[], double y[], size_t n);
  93.  
  94. xy_table
  95. make_xy_table (double x[], double y[], size_t n)
  96. {
  97.   xy_table t;
  98.   t.x = x;
  99.   t.y = y;
  100.   t.n = n;
  101.   return t;
  102. }
  103.  
  104. static int
  105. test_interp (
  106.   const xy_table * data_table,
  107.   const gsl_interp_type * T,
  108.   xy_table * test_table,
  109.   xy_table * test_d_table,
  110.   xy_table * test_i_table
  111.   )
  112. {
  113.   int status = 0;
  114.   size_t i;
  115.  
  116.   gsl_interp_accel *a = gsl_interp_accel_alloc ();
  117.   gsl_interp *interp = gsl_interp_alloc (T, data_table->n);
  118.  
  119.   gsl_interp_init (interp, data_table->x, data_table->y, data_table->n);
  120.  
  121.   for (i = 0; i < test_table->n; i++)
  122.     {
  123.       double x = test_table->x[i];
  124.       double y;
  125.       double deriv;
  126.       double integ;
  127.       double diff_y, diff_deriv, diff_integ;
  128.       gsl_interp_eval_e (interp, data_table->x, data_table->y, x, a, &y);
  129.       gsl_interp_eval_deriv_e (interp, data_table->x, data_table->y, x, a, &deriv);
  130.       gsl_interp_eval_integ_e (interp, data_table->x, data_table->y, 0.0, x, a, &integ);
  131.       diff_y = y - test_table->y[i];
  132.       diff_deriv = deriv - test_d_table->y[i];
  133.       diff_integ = integ - test_i_table->y[i];
  134.       if (fabs (diff_y) > 1.e-10 || fabs(diff_deriv) > 1.0e-10 || fabs(diff_integ) > 1.0e-10) {
  135.     status++;
  136.       }
  137.     }
  138.  
  139.   gsl_interp_accel_free (a);
  140.   gsl_interp_free (interp);
  141.  
  142.   return status;
  143. }
  144.  
  145. static int
  146. test_linear (void)
  147. {
  148.   int s;
  149.  
  150.   double data_x[4] = { 0.0, 1.0, 2.0, 3.0 };
  151.   double data_y[4] = { 0.0, 1.0, 2.0, 3.0 };
  152.   double test_x[6] = { 0.0, 0.5, 1.0, 1.5, 2.5, 3.0 };
  153.   double test_y[6] = { 0.0, 0.5, 1.0, 1.5, 2.5, 3.0 };
  154.   double test_dy[6] = { 1.0, 1.0, 1.0, 1.0, 1.0, 1.0 };
  155.   double test_iy[6] = { 0.0, 0.125, 0.5, 9.0/8.0, 25.0/8.0, 9.0/2.0 };
  156.  
  157.   xy_table data_table = make_xy_table(data_x, data_y, 4);
  158.   xy_table test_table = make_xy_table(test_x, test_y, 6);
  159.   xy_table test_d_table = make_xy_table(test_x, test_dy, 6);
  160.   xy_table test_i_table = make_xy_table(test_x, test_iy, 6);
  161.  
  162.   s = test_interp (&data_table, gsl_interp_linear, &test_table, &test_d_table, &test_i_table);
  163.   gsl_test (s, "linear interpolation");
  164.   return s;
  165. }
  166.  
  167. static int
  168. test_polynomial (void)
  169. {
  170.   int s;
  171.  
  172.   double data_x[4] = { 0.0, 1.0, 2.0, 3.0 };
  173.   double data_y[4] = { 0.0, 1.0, 2.0, 3.0 };
  174.   double test_x[6] = { 0.0, 0.5, 1.0, 1.5, 2.5, 3.0 };
  175.   double test_y[6] = { 0.0, 0.5, 1.0, 1.5, 2.5, 3.0 };
  176.   double test_dy[6] = { 1.0, 1.0, 1.0, 1.0, 1.0, 1.0 };
  177.   double test_iy[6] = { 0.0, 0.125, 0.5, 9.0/8.0, 25.0/8.0, 9.0/2.0 };
  178.  
  179.   xy_table data_table = make_xy_table(data_x, data_y, 4);
  180.   xy_table test_table = make_xy_table(test_x, test_y, 6);
  181.   xy_table test_d_table = make_xy_table(test_x, test_dy, 6);
  182.   xy_table test_i_table = make_xy_table(test_x, test_iy, 6);
  183.  
  184.   s = test_interp (&data_table, gsl_interp_polynomial, &test_table, &test_d_table, &test_i_table);
  185.   gsl_test (s, "polynomial interpolation");
  186.   return s;
  187. }
  188.  
  189.  
  190. static int
  191. test_cspline (void)
  192. {
  193.   int s;
  194.  
  195.   double data_x[3] = { 0.0, 1.0, 2.0 };
  196.   double data_y[3] = { 0.0, 1.0, 2.0 };
  197.   double test_x[4] = { 0.0, 0.5, 1.0, 2.0 };
  198.   double test_y[4] = { 0.0, 0.5, 1.0, 2.0 };
  199.   double test_dy[4] = { 1.0, 1.0, 1.0, 1.0 };
  200.   double test_iy[4] = { 0.0, 0.125, 0.5, 2.0 };
  201.  
  202.   xy_table data_table = make_xy_table(data_x, data_y, 3);
  203.   xy_table test_table = make_xy_table(test_x, test_y, 4);
  204.   xy_table test_d_table = make_xy_table(test_x, test_dy, 4);
  205.   xy_table test_i_table = make_xy_table(test_x, test_iy, 4);
  206.  
  207.   s = test_interp (&data_table, gsl_interp_cspline, &test_table, &test_d_table, &test_i_table);
  208.   gsl_test (s, "cspline interpolation");
  209.   return s;
  210. }
  211.  
  212.  
  213. static int
  214. test_akima (void)
  215. {
  216.   int s;
  217.  
  218.   double data_x[5] = { 0.0, 1.0, 2.0, 3.0, 4.0 };
  219.   double data_y[5] = { 0.0, 1.0, 2.0, 3.0, 4.0 };
  220.   double test_x[4] = { 0.0, 0.5, 1.0, 2.0 };
  221.   double test_y[4] = { 0.0, 0.5, 1.0, 2.0 };
  222.   double test_dy[4] = { 1.0, 1.0, 1.0, 1.0 };
  223.   double test_iy[4] = { 0.0, 0.125, 0.5, 2.0 };
  224.  
  225.   xy_table data_table = make_xy_table(data_x, data_y, 5);
  226.   xy_table test_table = make_xy_table(test_x, test_y, 4);
  227.   xy_table test_d_table = make_xy_table(test_x, test_dy, 4);
  228.   xy_table test_i_table = make_xy_table(test_x, test_iy, 4);
  229.  
  230.   s = test_interp (&data_table, gsl_interp_akima, &test_table, &test_d_table, &test_i_table);
  231.   gsl_test (s, "akima interpolation");
  232.   return s;
  233. }
  234.  
  235.  
  236. int 
  237. main (int argc, char **argv)
  238. {
  239.   int status = 0;
  240.  
  241.   argc = 0;    /* prevent warnings about unused parameters */
  242.   argv = 0;
  243.  
  244.   status += test_bsearch();
  245.   status += test_linear();
  246.   status += test_polynomial();
  247.   status += test_cspline();
  248.   status += test_akima();
  249.  
  250.   exit (gsl_test_summary());
  251. }
  252.